Theory of Diffusion Controlled Growth 
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We present a new theoretical framework for Diffusion Limited Aggregation and associated Dielec- 
tric Breakdown Models in two dimensions. Key steps are understanding how these models interrelate 
when the ultra-violet cut-off strategy is changed, the analogy with turbulence and the use of log- 
arithmic field variables. Within the simplest, Gaussian, truncation of mode-mode coupling, all 
properties can be calculated. The agreement with prior knowledge from simulations is encouraging, 
and a new superuniversality of the tip scaling exponent is both predicted and confirmed. 

PACS numbers: 61.43.Hv,47.53.+n 



Diffusion Limited Aggregation (DLA) has accumulated 
an enormous literature since Witten and Sander first in- 
troduced their simulation model of a rigid cluster growing 
by the accretion of dilute diffusing particles fjj]. The im- 
portance of the model is that it encompasses a range of 
problems where growth or interfacial advance is governed 
by a conserved gradient flux, that is the local interfacial 
velocity is given by 



v 2 ^ = 0, 



^interface 



(1) 



where for DLA 77 = 1 Q . The generalisation to a range of 
positive 77 was introduced by Niemeyer, Pietronero and 
Wiesmann || to model dielectric breakdown patterns, 
but in this letter we exploit it to support proposed equiv- 
alences between models with significantly different ultra- 
violet cut-off mechanism. Theoretical interest has been 
fuelled by the fractal and multifractal ||, || scaling prop- 
erties of the clusters produced, with controversial claims 
H |7], |f| (and counter-claims || [U^ |l]] ]) of anomalous 
scaling, and by the longstanding absence here resolved of 
an overall theortical framework to understand the prob- 
lem. 

The presence of a cut-off lcngthscale a below which the 
physics dictates smooth growth is a crucial ingredient of 
DLA; it is known that otherwise infinitely sharp cusps de- 
velop in the interface within finite time . In DLA this 
cutoff is fixed and set by the size of accreting particles, 
but there are other problems where it is set in a more sub- 
tle dynamical way by the surface boundary conditions on 
the diffusion field. In dendritic soldification this comes 
about through competition between surface energy and 
diffusion kinetics (with r\ = 1), leading to 



a cx \d n 4>[ 



(2) 



with m = 1/2 at least for those tips not in retreat p3[ . 
In terms of m, simple DLA corresponds to m = 0, and 
in the theory below in two dimensions we will map onto 
the case where a is such that each growing tip has fixed 
integrated flux, corresponding to m = 1. 

It is central to fractal (and multifractal) behaviour in 
DLA that the measure given by the diffusion flux [den- 
sity] d n (j) onto the interface has singularities |4|, |), such 



that the integrated flux onto the growth within distance 
r of a singular point is given by 



H{r) 



(3) 



Applying this phenomenology to the scaling around 
growing tips, we can establish an equivalence between 
models at different 77 and m by requiring that the relative 
advance rates of different growing tips are matched. Con- 
sider two competing tips labelled 1,2, for two growths 
with the same overall geometry but growing governed 
by parameters (77, m) and [rf ,ml) respectively. For 
tip 1 we will have tip radius a\ and flux density ji 
which arc matched between the two different models by 
j'xa'^/af = jia^ 1 / ? and similarly for tip 2, whilst 
the two tips are interrelated by a' x j'™ = a^j'-J 11 and sim- 
ilarly for the unprimed quantities. If we insist that their 
advance velocities are in the same ratio in both models 
this requires (ji/j^) 1 ' = Ui/j^ , which forces the pa- 
rameter relation 



1 + m(l +a-d) 1 + m'(l +a-d) 



Tj 



(4) 



For the two models to be equivalent in the relative ve- 
locities of all tips requires their parameters be related 
as above, where a is the singularity exponent associated 
with growing tips which we take to be the same as we 
are matching the geometry at scales above the cut-offs. 

Although we have not strictly proved the equivalence 
of the models related above, we have shown that any 
such relationship must follow Eq. (ji|) and we will assume 
in the rest of this letter that this equivalence holds. All 
such models are then classifiable in terms of a convenient 
reference such as 770, the equivalent 77 when m = 0, cor- 
responding to the original Dielectric Breakdown Model. 
For example dendritic solidification with 77 = 1 and 
rn = 1/2 corresponds to 770 = 3+ ^_ d ■ it is thus not 
equivalent to DLA, but to another member of the DBM 
class. Another puzzle resolved by our classification is 
a recent study showing conflicting scaling between DLA 
and different limits of a 'laplacian growth' model H] . In 
the present terminology the latter model corresponds to 
m = —1 and its two limits of low and high coverage of the 
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FIG. 1: The size fluctuations ("output noise") A out = 
(SN/N) 2 , measured at fixed radius R for DLA clusters grown 
with off- lattice noise reduction at various ("input") noise 
levels A[ n . For low A ln , A out self organizes from below in the 
manner of a turbulent system. 



growing surface per growth step have 77 = 3 and 77 = 1 re- 
spectively. Using a = 0.7 (see below) these map through 
Eq. (Q) into 770 = 2.31 and 770 = 0.77 respectively, so the 
way their scaling brackets that of DLA is quite expected. 

DLA and DBM have hitherto been regarded as models 
in statistical physics, in that the local advance rate in 
Eq. (Q) has been implemented as the probabilty per unit 
time for the growth locally to make some unit of advance, 
entailing an inherent shot noise. Here we argue that dif- 
fusion controlled growth is a problem of turbulence type, 
with noise self-organising from minimal input. The data 
in Fig. |l| show how the relative fluctuations can approach 
their limiting value from below as well as from above. 

The new ideas above, that we can balance changing 
the cut-off exponent m by adjustment of rj, and that 
noise can be left to self organise, are the key to a new 
theoretical formulation of the problem, at least in two 
dimensions of space to which we now specialise. In two 
dimensions the Laplace equation in (|l|) can be solved in 
terms of a conformal transformation between the physical 
plane of z = x + iy and the plane of complex potential 
lu = (p + i9 in which we take the growing interface to 
be mapped into the periodic interval 9 = [0, 2tt),<^ = 
and the region outside of the growth mapped onto <j> > 0. 
Then adapting reference |T^ |, we have for the dynamics 
of the interface following Eq. (|l|) , 



dz{6) 
dt 





de 


1+?) - 




dz 





(5) 



The linear operator V is most simply described 
in terms of Fourier transforms: T^Efe e ~ lke fk] — 



-ike 



fk 



fo + 2 I]f=i e tk6 fk, where we have 
introduced here an upper cut-off wavevector K . It is eas- 
ily shown that on scales of 9 greater than K~ x a smooth 
interface is linearly unstable with respect to corrugation 



for 77 > (the Mullins-Sekerka instability (L5)), whereas 
for scales of 9 less than K~ Y the equation drives smooth 
behaviour (corresponding locally to the case 77 = — 1). 
This cutoff on a scale of 9, the cumulative integral of 
flux, corresponds in terms of tip radii and flux densities 
to aj « K~ , that is an m — 1 cutoff law. Thus the pa- 
rameter 77 in Eq. (^|) is more specifically 771 = 0:770, using 
Eq. (g) with d = 2. 

We have made a numerical test of Eq. (JsJ) and the 
equivalence (g), with disorder supplied only through the 
initial condition, by applying them to the case of growth 
along a channel with periodic boundary conditions (cylin- 
der) . For this case analyticity of the conformal map re- 
quires that z(9) = i9 + J2k>o z k e ~ lk9 an d the overall ad- 
vance rate of the growth reduces to ^f- — ( 

which we can compare to the expected scaling of tip ve- 
locity with the cutoff, v ~ K^~ a ^^ a . It is convenient 
to change variables to ip = (dz/d9)~( 1+m ^ 2 , in terms of 
which we obtain 



d - 



(6) 



where y = (1 +771) /2 and the tri- linear form of the RHS 
enables us to compute numerically the motion within a 
purely Fourier representation. Figure ^| shows the mea- 
sured variation of J2j<k IV^'I 2 vs k m : this is expected to 
exhibit a power law with exponent (1/a — 1) and the 
observed slope plotted in this way is surprisingly inde- 
pendent of 771 . 

The most important result of our numerical study of 
Eq. is that this clearly does self-organise into sta- 
tistical scaling behaviour, given disorder from only the 
initial conditions. However the numerical results are also 
remarkable, as we obtain a « 0.74 ± 0.02 with no signif- 
icant dependence on 771 in the range studied. This not 
only agrees reasonably with the value a = D — 1 = 0.71 
known from large direct simulations of DLA |l6|, |l7| , but 
also appears to imply a deeper unversality which we will 
see is replicated in our analytic theory below. 

We now turn to a theoretical analysis of Eq. (j^), for 
which a primary requirement is that we must obtain re- 
sults explicitly independent of the cut-off as K — * 00. 
This is hard because we have already seen that the mean 
advance rate of the interface diverges as a power of K, 
and on fractal scaling grounds one would expect the same 
divergent factor to appear in the rate of change of simple 
variables such as z k or ip^. One can of course take ra- 
tios of rates of change and look to order terms such that 
divergences cancel, but to make this work we have been 
forced to introduce yet another change of variables, 



dz 

-i-QQ = exp(-A(0)) = cxp 




-ike 



(7) 



which corresponds to Fourier decomposing the logarithm 
of the flux density. The key to the success of these vari- 
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FIG. 2: Cumulative contribution to the mean growth velocity 
plotted against wavevector as fc' )1 with logarithmic scales. The 
data are (bottom to top) for rji — 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 and 
all exhibit a common power law slope 1/a-l* 0.35 ± 0.04 
per the guidelines shown. 



ables is that they decompose the flux density itself mul- 
tiplicatively and, as we shall see, quite naturally capture 
its multifractal behaviour. In terms of these 'logarithmic 
variables' the equation of motion becomes 



dt 



j<k 



P 



y(A+A) 



2k 



,t/(A+A) 



(8) 



where subscripts on bracketed expressions imply the tak- 
ing of a Fourier component, by analogy with A&. The 
advance rate of the mean interface is given in these vari- 
ables by = (e^ A+A )) . 

Now let us suppose some ignorance of the initial 
conditions and describe the system in terms of a joint 
probability distribution over the A&, and let us denote 
averages over this [unknown] distribution by (..). Wc 
can in principle determine the distribution through its 
moments, whose evolution we now compute. For sim- 
plicity we assume translational invariance with respect 
to 9, so that only moments of zero total wavevector need 
be considered, of which the lowest gives: ^ (AfeA^) = 



y(A+A)\ 



2k 



A fc ef+ X > V 



V j<fc x 

(c. conj.). All of the higher moments lead 

to the same form of averages on the RHS, 

^multinomial(A, \)e y ( x+ ^\ and all of these terms 
are conveniently expressed in terms of cumulants 
@, using the identities (Xe w ) / (e w ) = (Xe w ) , 
(XYe w )/(e w ) = (XYe w ) c + (Xe w ) 9 (Ye w ) c , etc. 
The key helpful feature is that the expressions we require 
all naturally divide by one factor of (^e y( - x+x ^ = ^ (zq), 
which is what we sought in order to remove divergences. 

To obtain tractable results we need to introduce some 
closure approximation(s) and we present here the sim- 
plest, neglecting all cumulants higher than the second, 



equivalent to assuming a joint Gaussian distribution (of 
zero mean) for A. This is entirely characterised by its sec- 
ond moments S(k) = (AfcA^) which by Eq. (H) we find 
evolve according to dS(k) / d (z ) — —kS(k) — y*kS(k) 2 — 
2y 2 J2j<k jSti)S(k) + 2ykS(k). This in turn approaches 
a stable steady state solution where 



S(k) 



2tf-l = 



k odd; S{k) = , k even. (9) 



The absence of even k is readily interpreted in terms of 
the dominance of one major finger and one major fjord. 

Within the Gaussian approximation and its predicted 
variances we can now compute all [static] proper- 
ties of diffusion controlled growth, in a channel and (see 
later discussion) also in radial geometry. The multifractal 
spectrum of the harmonic measure follows from comput- 
ing the general moment [|| /||f | T ^ = ^e (A+A)r/2 \ = 

exp(V 2 /4£f S(k)) 



q(r) = 1 + r + t 2 



leading to 



2(l + ryir 



(10) 



and it is easy to see that any closure scheme based on 
keeping cumulants of A up to some finite order leads 
to a polynomial truncation of q(r). From the Legen- 
dre Transform of the inverse function r(q) we obtain the 
corresponding spectrum of singularities, 



/(«) 



1 
















a 









(11) 



which in Fig. |3| is compared to measured data for DLA 
p9| , which later measurements |2(| reinforce. For the 
region of active growth a < 1 (q > 0) the theory is quan- 
titatively accurate. At a = 1 it conforms to Makarov's 
theore m plf , and in contrast to the Screened Growth 
Model p2| it does this without adjustment. For a > 1 
the spectrum only qualitatively the right shape, and for 
such screened regions our equations based on tip scaling 
may not hold. 

Although the multifractal moments depend signifi- 
cantly on the input parameter 771 , the tip scaling ex- 
ponent a turns out to be independent of this and in 
close agreement with our numerical results. Matching the 
expected scaling of the mean velocity (as used to mea- 
sure a above) to that of the multifractal moment with 
t = —(1 + 771) leads direcly to a = 2/3 independent of 
ifi . This is a remarkable success for the Gaussian Theory 
to predict this hitherto unexpected result so closely. 

The multifractal spectrum suggests that the Gaus- 
sian approximation is good in the growth zone, so we 
have computed the penetration depth as a further test. 
For growth in the channel we define relative penetra- 
tion depth 5 as the standard deviation of depth 3?(z) 
along the chanel, computed over the harmonic measure, 
divided by the width of the channel. This leads to 
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FIG. 3: Multifractal spectra from the Gaussian theory (a = 
2/3), compared to measured values for DLA E3|. Agreement 
is excellent for the active region r > 0, a < 1, and there are 
no adjustable parameters. 



■/2 = E 



fe>0 



/2 where the re- 



quired averages can all be computed in the approximation 
of Gaussian distributed A. Using r\\ = 2/3 correspond- 
ing to DLA this leads to Stheory = 0.13, compared to 
Sdla = 0.14 from direct simulations of DLA growth in 
a periodic channel |L7"| . 

All of the new theory is readily extended to growth 
from a point seed in radial geometry. The multifractal 
spectrum turns out to be unchanged, in accordance with 
expectations from universality The penetration depth 
relative to radius gives Stheory = 0.20 for radial DLA, 
compared to our recently published extropolation from 
simulations, Sdla = 0.13 @. 

For DLA and its associated Dielectric Breakdown 
Models we have shown a theoretical framework which 
is complete in the sense that essentially all measurable 
quantities can be calculated. For amplitude factors such 
as the relative penetration depth there is no theoreti- 
cal precedent. For the full spectrum of exponents the 
advance over the Screened Growth Model is the elimi- 
nation of fitting parameters. For the exponent evti P we 
have in the Gaussian approximation a striking new result 
that this is independent of 77, which begs direct confirma- 
tion by (expensive) particle-based simulations. However 
for DLA in particular we have not yet improved on the 
best theoretical value of atip, which remains 1/V2R* 0.71 
from the Cone Angle Approximation p3] . 

Within DLA and DBM we look forward to calculat- 
ing more properties such as the response to anisotropy, 
which is fairly readily incorporated into our equations of 
motion. A more challenging avenue is to improve on the 
Gaussian approximation which we have used to obtain 
explicit theoretical results. Truncating at a cumulant of 
higher order than the second is hard, and more seriously 
it does not correspond to a positive (semi-)definite prob- 
ability distribution. An alternative route of improvement 



which we are exploring is closure at the level of the full 
multifractal spectrum. 

There are possibilities for wider application of ideas in 
this letter, where we have formulated DLA and DBM as 
a turbulent dynamics governed by a complex scalar field 
in 1+1 dimensions. Decomposing this field multiplica- 
tively (through Fourier representation of its logarithm) 
was the crucial step to obtain renormalisable equations 
and theoretical access to the multifractal behaviour, even 
though other representations offered equations of motion 
(^) with weaker non-linearity. It is natural to specu- 
late whether the same strategy might apply to turbulent 
problems more widely, where the key issue appears to be 
identifying suitable fields to decompose multiplicatively 
which are of local physical significance, and subject to 
closed equations of motion. 

This research has been supported by a Marie Curie 
Fellowship of the EC programme "Improving Human Po- 
tential" under contract number HPMF-CT-2000-00800. 
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